late season snow

# Import SNOTEL data 
load("data/snotel_cont_data.RData")

# Import seasonal snow zone polygons
#SSZ <- st_read('data/MODIS_Snow_Zones/SSZ_0cc.shp')
PSZ <- st_read('data/MODIS_Snow_Zones/PSZ_0cc.shp') 

## Reading layer `PSZ_0cc' from data source 
##   `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\PSZ_0cc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357223 ymin: 1459534 xmax: -437223.4 ymax: 3149702
## Projected CRS: NAD_1983_Albers

eco_L3 <- st_read('data/eco_regions/us_eco_l3.shp')

## Reading layer `us_eco_l3' from data source 
##   `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\eco_regions\us_eco_l3.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1250 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2356069 ymin: 272048.5 xmax: 2258225 ymax: 3172577
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version

# Group sites by site_id
cont_snotel_site <- cont_snow_data %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise()

## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.

# transform files to CRS:4326
snotel_sites <- st_as_sf(cont_snotel_site, coords = c('longitude', 'latitude'), crs = 4326)
PSZ_4326 <- st_transform(PSZ, crs = 4326)
eco_L3_4326 <- st_transform(eco_L3, crs = 4326)

# clip SNOTEL Sites to PSZ
sf::sf_use_s2(FALSE)

## Spherical geometry (s2) switched off

PSZ_snotel_sites = st_intersection(PSZ_4326, snotel_sites)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

#Assign Eco Region Data
#PSZ_snotel_sites <- st_join(PSZ_snotel_sites, left = FALSE, eco_L3_4326[c("US_L3NAME", "US_L3CODE", "L3_KEY")])
snotel_sites_eco <- st_join(snotel_sites, left = FALSE, eco_L3_4326[c("US_L3NAME", "US_L3CODE", "L3_KEY")])

## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

#Create Map of SNOTEL Sites based on persistent snow zone (PSZ)
mapview(snotel_sites, col.regions = "red") +
  mapview(PSZ_snotel_sites, col.regions = "blue")

# Import data from .RData
#load("data/PSZ_cont_snotel.RData")
load("data/cont_snotel.RData")

# Count the number of the 
eco_sntl_sites <- cont_snotel %>% 
  group_by(US_L3NAME) %>% 
  summarize(eco_sntl_sites = n_distinct(site_id)) %>% 
  filter(eco_sntl_sites >= 5)

# All SNOTEL Sites, Calculate spring days with new SWE and new SWE depth 
# filter by date (March 1st - 60/61, March 15th - 74/75, April 1st - 91/92, April 15th - 105/106, August 1st - 243/244)
# Filter data by date from 2000 to present and from March 1st to August 1st
# group by year and site
start_year = 2000
# start_day = 60 #March 1
# cont_snotel_spring_mar1 <- cont_snotel %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
#   mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
#   filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
#                                                    yday >= start_day+1 & yday <= 243+1)) %>% 
#   select(-c("X", "network", "description", "start", "end"))
# 
# start_day = 74 #March 15
# cont_snotel_spring_mar15 <- cont_snotel %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
#   mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
#   filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
#                                                    yday >= start_day+1 & yday <= 243+1)) %>% 
#   select(-c("X", "network", "description", "start", "end"))

start_day = 91 #April 1
cont_snotel_spring_apr1 <- cont_snotel %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
  filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
                                                   yday >= start_day+1 & yday <= 243+1)) %>% 
  select(-c("X", "network", "description", "start", "end"))

# start_day = 105 #April 15
# cont_snotel_spring_apr15 <- cont_snotel %>%
#   mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
#   mutate(newSWE = ifelse(snow_water_equivalent>lag(snow_water_equivalent)+3, snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
#   filter(state != "SD", year >= start_year, ifelse(leap_year(year) == FALSE, yday >= start_day & yday <= 243, 
#                                                    yday >= start_day+1 & yday <= 243+1)) %>% 
#   select(-c("X", "network", "description", "start", "end"))

# # Determine number of days with increased SWE per spring at each site (YEARLY - March 1)
# cont_snotel_site_yr_mar1 <- cont_snotel_spring_mar1 %>%
#   group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
#   summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
#   merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
#   mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "3/1")
# 
# # Determine number of days with increased SWE per spring at each site (YEARLY - March 15)
# cont_snotel_site_yr_mar15 <- cont_snotel_spring_mar15 %>%
#   group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
#   summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
#   merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
#   mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "3/15")

# Determine number of days with increased SWE per spring at each site (YEARLY - April 1)
cont_snotel_site_yr_apr1 <- cont_snotel_spring_apr1 %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
  summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
  merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
  mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "4/1") 

## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude', 'elev', 'US_L3CODE', 'US_L3NAME'. You can override
## using the `.groups` argument.

# # Determine number of days with increased SWE per spring at each site (YEARLY - April 1)
# cont_snotel_site_yr_apr15 <- cont_snotel_spring_apr15 %>%
#   group_by(site_id, site_name, state, year, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY) %>%
#   summarise(newSWE_days = sum(newSWE>0, na.rm = TRUE), newSWE = sum(newSWE, na.rm = TRUE)) %>% 
#   merge(., eco_sntl_sites, by = c("US_L3NAME")) %>% 
#   mutate(eco_sntl = paste0(L3_KEY, ' (N = ', eco_sntl_sites, ')'), start = "4/15")

# Merge the DF's into one large datafame for plotting
# cont_snotel_site_yr <- rbind(cont_snotel_site_yr_mar1, cont_snotel_site_yr_mar15, cont_snotel_site_yr_apr1, cont_snotel_site_yr_apr15)

# Determine average number of days with increased SWE per spring at each site in the PSZ (SITE - April 1)
cont_snotel_site_apr1 <- cont_snotel_site_yr_apr1 %>%
  group_by(site_id, site_name, state, latitude, longitude, elev, US_L3CODE, US_L3NAME, L3_KEY, eco_sntl_sites, eco_sntl) %>%
  summarise(mean_spring_days = mean(newSWE_days), med_spring_days = median(newSWE_days),
            mean_newSWE = mean(newSWE), med_newSWE = median(newSWE))

## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude', 'elev', 'US_L3CODE', 'US_L3NAME', 'L3_KEY',
## 'eco_sntl_sites'. You can override using the `.groups` argument.

# # View data
# mapview(PSZ_cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "med_spring_days",
#         layer.name = "AVG Sping Days (PSZ Sites)" , crs = 4326, grid = FALSE) +
#   mapview(eco_L3_4326, zcol = "L3_KEY")

# Spring SWE increases (After April 1)
spring_snow_days <- ggplot(cont_snotel_site_yr_apr1, aes(x = reorder(US_L3CODE,newSWE_days,na.rm = TRUE), y= newSWE_days, 
                                                         fill = reorder(eco_sntl,newSWE_days,na.rm = TRUE))) +
  geom_boxplot() +
  labs(x = "EPA Level III Ecoregion", y = "Days of Increased SWE", fill = "") +
  scale_y_continuous(breaks = seq(0,50,10), limits = c(0,50), expand = c(0.01, 0.01)) +
  theme_bw()
#spring_snow_days

ggplotly(spring_snow_days)

#ggplot_build(spring_snow_days)$data

# Elevation vs new SWE count (March 1)
elev_newsnow <- ggplot(cont_snotel_site_apr1, aes(x = med_spring_days, y= elev, 
                                                  color = reorder(eco_sntl,med_spring_days,na.rm = TRUE))) +
  geom_point() +
  labs(x = "Median Days of Increased SWE", y = "Elevation (m)", color = "") +
  scale_x_continuous(breaks = seq(0,30,10), limits = c(0,30), expand = c(0.01, 0.01)) +
  scale_y_continuous(breaks = seq(500,4000,1000), limits = c(500,3750), expand = c(0.01, 0.01)) +
  scale_color_manual(values=c("#F8766D", "#E58700", "#C99800", "#A3A500", "#6BB100",
                              "#00BA38", "#00BF7D", "#00C0AF", "#00BCD8", "#00B0F6",
                              "#619CFF", "#B983FF", "#E76BF3", "#FF67A4", "#FD61D1")) +
  theme_bw()
#elev_newsnow
ggplotly(elev_newsnow)

plot <- plot_grid(
  spring_snow_days + theme(legend.position="none"),
  elev_newsnow + theme(legend.position="none"),
  align = 'vh',
  labels = c("a.", "b."),
  hjust = -1,
  nrow = 1)

plot

legend <- get_legend(spring_snow_days + guides(fill = guide_legend(nrow = 5)) + 
               theme(legend.direction = "horizontal",
                     legend.justification="center",
                     legend.box.just = "bottom"))

plot_grid(plot, legend, ncol = 1, rel_heights = c(1, .3))

# # Elevation vs new SWE depth
# elev_newsnow <- ggplot(cont_snotel_site_apr1, aes(x = med_newSWE, y= elev, color = eco_sntl)) +
#   geom_point() +
#   labs(title = "Average Increase in SWE After April 1st (2000-2022)", x = "Median Depth of Increased SWE (mm)", y = "Elevation (m)", color = "Ecoregions (Number of SNOTEL Sites)")
# ggplotly(elev_newsnow)

# elev_newsnow <- plot_ly(cont_snotel_site, x = ~med_spring_days, y = ~elev, color = ~US_L3NAME, type = 'scatter', mode = 'markers', 
#                         hoverinfo = 'text',
#                         text = ~paste('</br> median new days: ', med_spring_days,
#                                       '</br> elev: ', elev,
#                                       '</br> site_name: ', site_name,
#                                       '</br> ecoregion: ', US_L3NAME,
#                                       '</br> state: ', state)) %>% 
#   layout(title = "Average Increased SWE Days After April 1st by Elevation (2000-2022)")
#  
# elev_newsnow

# elev_newsnow_depth <- plot_ly(PSZ_cont_snotel_site, x = ~med_newSWE, y = ~elev, color = ~US_L3NAME, type = 'scatter', mode = 'markers', 
#                         hoverinfo = 'text',
#                         text = ~paste('</br> median new SWE (mm): ', med_newSWE,
#                                       '</br> elev: ', elev,
#                                       '</br> site_name: ', site_name,
#                                       '</br> ecoregion: ', US_L3NAME,
#                                       '</br> state: ', state)) %>% 
#   layout(title = "Average Increased SWE Depth (mm) After April 1st by Elevation (2000-2022)")
#  
# elev_newsnow_depth